IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 Image Deblurring Using Derivative Compressed Sensing for Optical Imaging Application Mohammad Rostami, Student Member, IEEE, Oleg Michailovich, Member, IEEE, and Zhou Wang, Member, IEEE o 00 < O > o o> in l> o X Abstract — Reconstruction of multidimensional signals from the samples of their partial derivatives is known to be a standard problem in inverse theory. Such and similar problems routinely arise in numerous areas of applied sciences, including optical imaging, laser interferometry, computer vision, remote sensing and control. Though being ill-posed in nature, the above problem can be solved in a unique and stable manner, provided proper regularization and relevant boundary conditions. In this paper, however, a more challenging setup is addressed, in which one has to recover an image of interest from its noisy and blurry version, while the only information available about the imaging system at hand is the amplitude of the generalized pupil function (GPF) along with partial observations of the gradient of GPF's phase. In this case, the phase-related information is collected using a simplified version of the Shack-Hartmann interferometer, followed by recovering the entire phase by means of derivative compressed sensing. Subsequently, the estimated phase can be combined with the amplitude of the GPF to produce an estimate of the point spread function (PSF), whose knowledge is essential for subsequent image deconvolution. In summary, the principal contribution of this work is twofold. First, we demonstrate how to simplify the construction of the Shack-Hartmann interferometer so as to make it less expensive and hence more accessible. Second, it is shown by means of numerical experiments that the above simplification and its associated solution scheme produce image reconstructions of the quality comparable to those obtained using dense sampling of the GPF phase. Index Terms — Deconvolution, inverse problem, derivative com- pressive sampling, and Shack-Hartmann interferometer. I. Introduction THE necessity to recover digital images from their dis- torted and noisy observations arises in a multitude of practical applications, with some specific examples including image denoising, super-resolution, image restoration, and wa- termarking, just to name a few |[l)-|j4J. In such cases, it is standard to assume that the observed image v is formed by subjecting the original image u to convolution with a point spread functiorT] (PSF) i, followed by contamination by white Gaussian noise (WGN) v. Thus, formally. i-^ u + V. (1) While u and v can be regarded as general members of the signal space L2(ri) of real-valued functions on 51 C M?, the PSF i is normally a much smoother function, with effectively band-limited spectrum. As a result, the convolution with i has a destructive effect on the informational content of u, in which All the authors are with the Department of Electrical and Computer Engineering, University of Waterloo, N2L 3G1 Ontario. ' Note that, in optical imaging, this function is also referred to as an impulse transfer function |5 1. case V typically has a substantially reduced set of features with respect to u. This makes the problem of reconstruction of u from V a problem of significant practical importance 161. Reconstruction of the original image u from v can be carried out within the framework of image deconvolution, which is a specific instance of a more general class of inverse problems JT). Most of such methods are Bayesian in nature, in which case the information lost in the process of convolution with i is recovered by requiring the optimal solution to reside within a predefined functional class |8 1, |9|. Thus, for example, in the case when u is known to be an image of bounded variation, the above regularization leads to the famous Rudin-Osher-Fatemi reconstruction scheme, in which u is estimated as a solution to the following problem 1 10 1, 1 1 1 1 arg mm 2" \l+a |Vu| dxdy > (2) where a > is the regularization parameter. It should be noted that, if the PSF obeys J i dxdy ^ 0, the problem O) is strictly convex and therefore admits a unique minimizer, which can be computed by a spectrum of available algorithms pO) , pTj . A particularly non-trivial version of deconvolution is com- monly referred to as blind. In this case, the original image u is to be estimated without the knowledge of the PSF I?). In this paper, however, we follow the philosophy of hybrid deconvolution |12|, which takes advantage of any partial information on the PSF to improve the image reconstruction. Thus, in the algorithm described in this paper, the original image u will be recovered from v and some partial information on i. Optical (and, in particular, turbulent) imaging is unarguably the field of applied sciences from which the notion of decon- volution has originally emanated p3j-p7j. In short-exposure imaging, however, computational methods of image restoration are still superseded by adaptive optics. As recently as a decade ago, the use of adaptive optics would have been considered as the only practical option. Nowadays, however, with the advent of distributed cluster computing and GPU- based image processing, it seems to be time to revisit the cost- to-performance characteristics of the existing tools of adaptive optics. Thus, in this work, our focus is on a specific tool of adaptive optics, known as the Shack-Hartmann interferometer 1 16 1, p7| . Instead of completely excluding the interferometer from our measurement system, we propose to modify its construction through reducing the number of its local wave- front lenses. Although the advantages of such a simplification are immediate to see, its main shortcoming is obvious as well: the smaller the number of lenses is, the stronger is IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 the effect of undersampling and aliasing. Accordingly, to overcome this problem, we propose to augment the modified Shack-Hartmann interferometer by subjecting its output to the derivative compressed sensing (DCS) algorithm of 1 18 1. As it will be shown later in the paper, the PSF i is determined by a generalized pupil function P, which can be expressed in a polar form as P = A e^"^. While the amplitude A can be mea- sured via calibration or computed as a function of the aperture geometry, the phase is often influenced by environmental effects and hence it needs to be recovered from observations. It will be shown below that DCS is particularly well suited for reconstruction of </) from incomplete measurements of its partial differences. Such an estimate can be subsequently combined with A to yield an estimate of the PSF i, which can in turn be used by a deconvolution algorithm. Thus, the proposed method for estimation of the PSF and subsequent deconvolution of u can be regarded as a hybrid deconvolution technique, which comes to simplify the design and complexity of adaptive optics on the one hand, and to make the process of reconstruction of optical images as automatic as possible, on the other hand. The rest of the paper is organized as follows. Section II summarizes basic technical preliminaries. In Section III, we describe the SH interferometer as well as phase measurements in optical imaging. In Section IV, we explain DCS and our new approach to solve it. In Section V we describe deconvolution process to recover the original image. Experimental results are presented in Section VI, while Section VII finahzes the paper with a discussion and main conclusions. II. Technical Preliminaries In short exposure imaging, due to aberrations in the imaging system induced by, e.g., atmospheric turbulence, the impulse response of an optical imaging system is often unknown p9) . In order to better understand the setup under consideration, we first note that, in optical imaging, the PSF i is obtained from an amplitude spread function (ASF) h as i :— |ft,p. The ASF, in turn, is defined in terms of the generalized pupil function (GPF) P(x, y) as given by 1,20} K^.ii)^ 1 A. P{x,y)e rix^+yri) dxdy, (3) where Zi is the focal distance and A^, is the optical wavelength. Being a complex-valued quantity, P{x, y) can be represented in terms of its amplitude A{x, y) and phase (/)(x, y) as P{x,y):^A{x,y)e^'f''-'-''y\ (4) Here, the GPF amplitude A{x,y) (which is sometimes simply referred to as the aperture function) is normally a function of the aperture geometry. Thus, for instance, in the case of a circular aperture, A{x, y) can be defined as ||19| A{r) = ifr<f otherwise (5) where D denotes the pupil diameter. Thus, given (j){x,y), one could determine h and therefore i. Unfortunately, the phase (j){x, y) does not have an analytic expression, and it has to be measured in practice using such tools as the Shack-Hartmann interferometer (SHI) | fT6l . As will be discussed later in the paper, the SHI is capable of sensing the partial derivatives of (t>{x,y). Needless to say, in order to minimize the effect of aliasing on the estimation result, an accurate reconstruction of 0(a:, y) requires taking a fairly large number of the samples of V0(x, y) |21 1. In some applications, the number of sampling points (as defined by the number of local wavefront lenses) reaches as many as a few thousands. It goes without saying that reducing the number of lenses would have a positive impact on the SHI in terms of its cost and approachability. Alas, such a reduction is impossible without undersampling, which tends to have formidable effect on the overall quality of phase estimation. In this paper, to minimize the effect of undersampling, we exploit DCS |18|. As opposed to the classical compressed sensing (CCS) |22| , in addition to the sparsifing constraints, DCS also uses constraints which are intrinsic in the definition of partial derivatives. Using these additional constraints - which are called the cross-derivative constraints - allows sub- stantially improving the quality of reconstruction of (t)[x,y), as compared to the case of CCS-based estimation. III. Shack-Hartmann Interferometer (SHI) As it was mentioned earlier, the SHI is typically used to measure the gradient V0(x, y) of the GPF phase 0(a;, y), from which the values of the latter can be subsequently estimated. To this end, the unknown phase (f){x, y) is assumed to be expandable in terms of some basis functions {Zk}'^^^, viz. ^ ' oo (l^{x,y)^^akZk{x,y), (6) fc=0 where the representation coefficients {ak}kLo are assumed to be unique and stably computable. Note that, in this case, the datum of {flfej^Q uniquely identifies (f){x,y), while the coefficients {a/d^g can be estimated due to the linearity of (|6]l which suggests V(?!)(x, y) = ^ak ^Zk{x, y), (7) fe=0 The most frequent choice of {^fe}^o ^'^ adaptive optics is the Zernike polynomials (aka Zemike functions) [20 1. These polynomials constitute an orthonormal basis in the space of square-integrable functions defined over the unit disk in M^. Zernike polynomials can be subdivided in two subsets of the even Z™ and odd Z,7™ Zernike polynomials which have very convenient analytical definitions as given by Z::^ip,^)^R^{p)cos{mip) (8) (9) where m and n are nonnegative integers with n > m, < If < 271 is the azimuthal angle, and < p < 1 is the radial distance. The radial polynomials _R™ are defined as (n-m)/2 Kip) (-!)'=(« -fc)! ^ A:!((n + m)/2-/s)!((n-m)/2-fc)! 11-2 k (10) IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 Fig. 1. An example of a 10 X 10 SHI array on a circular aperture. The shading indicates those blocks (i.e., lenses) which are rendered idle. Note that, since the Zernike polynomials above are defined using polar coordinates, it makes sense to re-express the phase (p and its gradient in the polar coordinate system as well. (Technically, this would amount to replacing x and y in (|6]l-(j7]i by p and ip, respectively.) Moreover, due to the property of the Zernike polynomials to be an orthonormal basis, the representation coefficients {a^l^Q in in ((6])-(|7]i can be computed by orthogonal projection, namely a-k (t){p,(p) Zk{p,(p) pdpdip (11) In practice, however, (j){p, (p) is unknown and therefore the coefficients {flfej^o n^ed to be estimated by other means. Thus, in the case of the SHI, the coefficients can be estimated from a finite set of discrete measurements of V0(p, (p). The main function of the SHI is to acquire discrete measure- ments of V0 by means of linearization. The linearization takes advantage of subdividing a (circular) aperture into rectangular blocks with their sides formed by a uniform rectangular lattice. An example of such a subdivision is shown in Fig. [Tlfor the case of a 10 X 10 lattice grid. Subsequently, it is assumed that the grid is sufficiently fine to approximate a restriction of the phase <j) to each of the above blocks by a linear function. This results in a piecewise linear approximation of 0, whose accuracy improves asymptotically when the lattice size goes to infinit>|^ Formally, let n := {{x,y) e M^ | x^ + y^ < jj2y be a circular aperture of radius D and 5 = {(x, y) e M^ | max{|a;|, \y\} < D} be a square subset of R-^ such that fi C 5. Then, for each polar coordinate (p, ip) £ ft and an N x N grid of square blocks of size 2D/N x 2D/N, the phase (j) can be expressed as 4>(x,y) ~ ax + by + c^ (12) for all (x, y) in a neighbourhood of (p cos (p,p sin ip). The approximation in ([T2| suggests that V0(x,y)« [a, 6]^ (13) ^More rigorously, one can show that, as long as </> is uniformly continuous, its piecewise linear approximation converges uniformly, as the grid size goes to infinity. LoisAnij C CD Clip Shack-Haattraan Pattern Fig. 2. Basic structure of the SHI and a resulting pattern of the focal points. where (•)-^ denotes matrix transposition. While c in ( 12 1 can be derived from boundary conditions, coefficients a and b should be determined via direct measurements. To this end, the SHI is endowed with an array of small focusing lenses, which are supported over each of the square blocks of the discrete grid. In the absence of phase aberrations, the focal points of the lenses are spatially identified and registered using a high-resolution CCD detector, whose imaging plane is aligned with the plane of the focal points. Then, when the wavefront gets distorted as a result of, e.g., atmospheric turbulence, the focal points are "pushed" towards new spatial positions, which can also be pinpointed by the same detector. The resulting displacements can therefore be measured and subsequently related to the values of V0 at corresponding points. To explain how the above procedure can be performed, additional notations are in order. Let ft^ denote a finite set of spatial coordinates defined according to n. Xd yd \{xd,yd) e n -D+^(i + ^], t^O,l,...,N~l (14) and xl + yj<D^y Note that the set fid can be thought of as a set of the spatial coordinates of the geometric centres of the wavefront lenses, restricted to the domain of aperture O. Under some reasonable assumptions, one can then show 1 ,23) that the focal displacement A{x,y) — [Ax{x, y), Ay{x,y)Y' measured at some {x, y) G fid is related to the value of V(/)(x, y) as given by V0(a;, y) « -A^{x, y), V(a;, y) G ^d, (15) where / is the focal distance of the wavefront lenses. Such a measurement setup is depicted in Fig.|2]along with an example of measured focal points. Now, given a total of M :— #51^ measurements of V0 over rijj, one can try to recover a useful approximation of (j) over the whole Vl in the form of a projection of 4> on the linear subspace spanned by all Zernike polynomials up to the IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 order L inclusive. In this case, it is possible to estimate the representation coefficients {ak\k=a ^^^ solution of conditions of compressed sensing far superseded those of the subject to appropriate boundary conditions. It is worthwhile noting that (16 1 can be rewritten in a vector-matrix form as IZa-dl 2 2j s.t. a h 0, (17) where Z is a 2Af x i + 1 matrix composed of the values of the partial derivatives of the Zernike polynomials, d is a measurement (column) vector of length 2M, and a — [oq, fli, . . . , ai]^ is a vector of the representation coefficients of (j). The constraint a )^ in (17 1 is supposed to further regularize the solution by forcing a to be a member of some convex cone as well. Thus, for example, if the mean value of (j) can be assumed to be equal to zero, the solution to ( 17 1 can be computed as a=Z#d, (18) where Z* denotes the pseudo-inverse of Z, whose definition is unique and stable as long as the row-rank of Z is greater or equal to L + 1 (hence suggesting that 2M > L + 1). Having estimated a, the phase can be approximated as <P{P,'~P) ~ ^akZk(p,(p). (19) k=0 A higher accuracy of phase estimation requires using higher- order Zernike polynomials, which in turn necessitates a pro- portional increase in the number of wavefront lenses. More- over, as required by the linearization procedure in the SHI, the lenses have to be of a relatively small sizes (sometimes, on the order of a few microns), which may lead to the use of a few thousand lenses per one interferometer. Accordingly, to simplify the construction and to reduce the cost of SHIs, we propose to substantially reduce the number of wavefront lenses, while compensating for the induced information loss through the use of DCS, which is detailed next. IV. Derivative Compressive Sampling A. Classical Compressed Sensing Central to signal processing is the Shannon-Nyquist theorem p4) , which specifies conditions on which a band-limited signal can be stably and uniquely recovered from its discrete measurements. However, in around 2005, a different sampling theorem was formed which, in some cases, abrogates the fundamentals of its predecessor. This new theory, nowadays known as compressed sensing (aka compressive sampling), asserts that signals, which admit a sparse representation in a predefined basis/frame, can be recovered from their discrete measurements, whose number is proportional to the ^o-norm of the coefficients of the sparse representation. In such a case, the sparser the representation of the signal is, the smaller can be the number of measurements required for signal reconstruction. As a result, cases are numerous in which the Shannon-Nyquist sampling |22|, |25|. Despite its widespread success in countless applications, the theory of compressed sensing is not entirely free of limitations. One of such limitations stems from the necessity to use a non- linear decoder. Indeed, while in the case of classical sampling, the reconstruction of a time-domain signal is implemented through linear interpolation (i.e., linear filtering), in the case of compressed sensing, the reconstruction involves solution of a convex optimization problem. Specifically, let us consider a typical setup of classical compressed sensing (CCS), in which y £ M'" represents an observed version of x e M", related according to y = *x, (20) where ^ G M™^" is an observation (sampling) matrix with n > m. The recovery of x from y based on ( [20| is impossible to implement in a unique and stable way, unless it is known that X is sparse and hence has a relatively low value of ||x||o. In such a case, if the sampling matrix ^ satisfies the restricted isometry property (RIP) f22|, p5) with respect to a certain class of sparse signals to which x is assumed to belong, then CCS recovers x as a solution to | |26l , | |27l X = argmin{||x'||i | ^x' = y} , (21) x' which is a convex minimization problem, which is straightfor- ward to reformulate in terms of linear programming. Moreover, in the case when the measurements y are error-prone, a more robust version of CCS is to recover x as given by X — arg mm {llx'l |*x' y\\l <^}. (22) where e > is a parameter controlling the size of the noise. Moreover, it was shown in 1221, ||25l, that the estimation error in the signal reconstructed according to ( 22 1 can be bounded by a linear function of e. This implies robustness of the CCS reconstruction towards the presence of measurement noise. It should be finally noted that the optimization problem (22 1 can be reformulated in its alternative Lagrangian form. in which case one can find 1 X — arg mm *x' y||^ Allx'l (23) where A > is an optimal Lagrange multiplier p6) . In what follows, it is assumed that an optimal value of A is known. (For more details on this subject, the reader is referred to p6| as well as to the later sections of this paper). B. Derivative Compressed Sensing (DCS) Let the partial derivatives of evaluated (by means of the SHI) at the points of set ilij be column-stacked into vectors fj; and iy of length M = i^fld- In what follows, the partial derivatives i^ and fj, are assumed to be sparsely representable by an orthonormal basis in M*^. Representing such a basis by an M X M unitary matrix W, the above assumption suggests the existence of two sparse vectors c^, and Cy such that Wcx and Wcy provide accurate approximations of the partial derivatives of the original phase (j) evaluated over fid- IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 Now, the simplification of the SHI proposed in the current paper amounts to reducing the number of wavefront lenses to a minimum. Formally, such a reduction can be described by two n X M sub-sampling matrices '^^ and '^y, where n < M. Specifically, let h^ '■— 'i'xix and bj, :— 'i'yiy be incomplete (partial) observations of i^ and fj^, respectively. Then, the noise-free counterparts of the partial derivatives can be approximated by Wc^ and Wc* respectively, where < = argmin |i||*,W^< - b,||2 + A,||<||i| (24) where p^*) is a vector of Bregman variables (or, equivalently, augmented Lagrange multipliers) and (5 > is a user-defined parameteipl and argmin<; -W'i'yWc'y - hy\\l ^v\W% Vll^yl (25) for some Aa;,Ay > 0. Moreover, in the case when A^^ — Xy, the above estimates can be combined together To this end, let c = [Cx,Cyf, b = [K,hyf, and A = diag{*^M^, *yTy} e M2"x2M^ Then, argmm ■^ -||Ac' b||2 + A||c'||i (26) The c-update step in ( 32 1 has the format of a standard basis pursuit de-noising (BPDN) problem f29l, which can be solved by a variety of optimization methods [301. In the present paper, we used the FISTA algorithm of |31) due to the simplicity of its implementation as well as for its remarkable convergence properties. It should be noted that the algorithm does not require explicitly defining the matrices A and B. Only the operations of multiplication by these matrices and their transposes need to be known, which can be implemented in an implicit and computationally efficient manner Once an optimal c* is recovered, it can be used to estimate the noise-free versions of i^ and fy as WT^c* and WTyC*, respectively. These estimates can be subsequently passed on to the fitting procedure of Section III to recover the values of (j), which, in combination with a known aperture function A, provide an estimate of the PSF i as an inverse discrete Fourier transform of the autocorrelation of P = AeJ'^. Algorithm 1 below summarizes our method of estimation of the PSF. where A = A, = Xy. In this form, the problem (|26f is identical Algorithm 1: PSF estimation via DCS to (23 I and hence it can be solved by a variety of optimization algorithms p6), p7). The DCS algorithm augments CCS by subjecting the mini- mization in (|26| to an additional constraint which stems from the fact that |18| dx dy dy dx ' which is valid for any two times continuously differentiable 4i{x^ y). Thus, in particular, the constraint implies the existence of two partial differences matrices D^ and Dy which obey DJy = Dyix. (28) Consequently, if T^. and Ty are the matrices satisfying T^c — 1) Data: h^, by, and A > 2) Initialization: For a given transform matrix W and matrices/operators '^x, ^y, D^, Dy, T^ and Ty, preset the procedures of multiplication by A, A^ , B and B^ . 3) Phase recovery: Starting with an arbitrary c^^^ and p'-''-' = 0, iterate ( 32 1 until convergence to result in an Cx and TyC = Cy, respectively, then (31 1 can be re-expressed in terms of c as (29) or ^y-^x^ — J-^x-^y^ Bc = 0, (30) where B :— DyTx — DxTy. Thus, DCS solves the constrained minimization problem as given by c* = argmin |i||Ac' ~h\\l + A||c'||i| , (31) s.t. Be' = optimal c*. Use the estimated (full) partial derivatives WTxC* and WTyC* to recover the values of (p over H.. 4) PSF estimation: Using a known aperture function A, compute the inverse Fourier transform of P = A eJ^ to result in a corresponding ASF h. Estimate the PSF i as ^^\h\\ The estimated PSF can be used to recover the original image u from V through the process of deconvolution as explained in the section that follows. V. Deconvolution The acquisition model ([TJ can be rewritten in an equivalent operator form as given by V ~ H{u} + V, (33) A solution to ( 3 1 1 can be found, for instance, by means of the Bregman algorithm L28J, in which case c* is approximated by a stationary point of the sequence of iterations produced by ^(t+i) ^ arg miiic' (211^'- - 1^112- \Ac'-h\\ +A||c'||i I Be' .(*) IB} (32) where Ti, denote the operator of convolution with the estimated PSF i. Note that, in this case, the noise term v accounts for both measurement noise as well as the inaccuracies related to estimation error in i. The deconvolution problem of finding a useful approxima- tion of u given its distorted measurement v can be addressed in many way, using a multitude of different techniques pTj-p3). In this work, we use the ROF model and recover a regularized ^In this work, we use 5 = 0.5. IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 approximation as a solution of u* = argniin <^ -\\n{u} - v\\l + 7 ||u|JTy u I Z (34) where ||u||tv = J J\Wu\dxdy denotes the total variation (TV) semi-norm of u. One computationally efficient way to solve ( (34] i is to substitute a direct minimization of the cost function in (|34]) by recursively minimizing a sequence of its local quadratic majorizers |31|. In this case, the optimal solution u* can be approximated by the stationary point of a sequence of intermediate solutions produced by w («) ,(*) M'H*{w--H{w(*'}} ,(t+i) ^ argmin„ {^\\u ~ u.(*)||2 + 7 \\u\\tv} , (35) where H* is the adjoint of H and fi is chosen to satisfy /i > IIH*?^!!. In this paper, the TV denoising at the second step of ([35| has been performed using the fixed-point algorithm of Chambolle fTTl. The convergence of ( |35] l can be further improved by using the same FISTA algorithm of pT| . The resulting procedure is summarized below in Algorithm 2. Algorithm 2: TV deconvolution using FISTA 1) Initialize: Select an initial value m'^^; set y'"^ = u^^^^ and t(o) = 1 2) Repeat until convergence: w (t+i) y^*'> + ^lH* {v - H{y^*^}} argmin„{i||M-u;(*)|||+7||u|JTy} (*+i) = 0.5 (1 + v'1 + 4(t(*))2) ,(*+!) In summary. Algorithms 1 and 2 represent the essence of the proposed algorithm for hybrid deconvolution of optical images. The next section provides experimental results which further support the value and applicability of the proposed methodology. VI. Results To demonstrate the viability of the proposed approach, its performance has been compared against reference methods. The first reference method used a dense sampling of the phase (as it would have been the case with a regular SHI), thereby eliminating the need for a CS-based phase reconstruction. The resulting method is referred below to as the dense sampling (DS) approach. Second, to assess the importance of incorpora- tion of the cross-derivative constraints, we have used both CCS and DCS for phase recovery. In what follows, comparative results for phase estimation and subsequent deconvolution are provided for all the above methods. A. Phase recovery To assess the performance of the proposed and reference methods under controllable conditions, simulation data was used. The random nature of atmospheric turbulence necessi- tated the use of statistical methods to model its effect on a (a) (b) (c) Fig. 3. An example of a simulated phase (j> (subplot (a)) along with its partial derivatives w.r.t. x (subplot (b)) and y (subplot (c)). wavefront propagation. Specifically, in this paper, the effect of atmospheric turbulence has been described using the modified Von Karman PSF model |34|. A typical example of a GPF phase (p is shown in subplot (a) of Fig. [5] In the shown case, the size of the phase screen was set to be equal to 10 x 10 cm, while the sampling was performed over a 128 x 128 uniform lattice (which would have corresponded to the use of 16384 lenses of a SHI). The partial derivatives d(j)/dx and d(j)/dy are shown in subplots (b) and (c) of Fig. [3] respectively. In the present paper, the subsampling matrices '^^ and ^y were obtained from an identity matrix / through a random subsampling of its rows to result in a required compression ratio r (to be specified below). To sparsely represent the partial derivatives of (j)^ ^ was defined to correspond to a four- level orthogonal wavelet transform using the nearly symmetric wavelets of Daubechies with five vanishing moments |35J. To demonstrate the value of using the cross derivative con- straint for phase reconstruction, the CCS and DCS algorithms have been compared in terms of the mean square errors (MSE) of their corresponding phase estimates. The results of this comparison are summarized in Fig. |4]for different compression ratios (or, equivalently, (sub)sampling densities) and SNR = 40 dB. As expected, one can see that DCS results in lower values of MSE as compared to CCS, which implies a higher accuracy of phase reconstruction. Moreover, the difference in the per- formances of CCS and DCS appears to be more significant for lower sampling rates, while both algorithms tend to perform similarly when the sampling density approaches the DS case. Specifically, when the sampling density is equal to r = 0.3, DCS results in a ten times smaller MSE as compared to the case of CCS, whereas both algorithms have comparable performance for r — 0.83. This result characterizes DCS as IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 Q H n — [] — ei Fig. 4. MSB of phase reconstruction obtained with different methods as a function of r. Here, the dashed and solid lines correspond to CCS and DCS, respectively, and SNR is equal to 40 dB. a better performer than CCS in the case of relatively low (sub)sampling rates. Some typical phase reconstruction results are shown in Fig.|5] whose left and right subplots depict the phase estimates obtained using the CCS and DCS algorithms, respectively, for the case of r = 0.5. To visualize the differences more clearly, error maps for both methods are shown in subplot (c) and (d)(note for better visualization error map values are multiplied by factor of 10.). A close comparison with the original phase (as shown in subplot (a) of Fig. [3]) reveals that DCS provides a more accurate recovery of the original (f), which further supports the value of using the cross-derivative constraints. In fact, exploiting these constraints effectively amounts to using additional measurements, which are ignored in the case of CCS. To investigate the robustness of the compared algorithms towards the influence of additive noises, their performances have been compared for a range of SNR values. The results of this comparison are summarized in Fig. [6] Since the cross- derivative constraints exploited by DCS effectively restrict the feasibility region for an optimal solution, the algorithm exhibits a substantially better robustness to the additive noise as compared to the case of CCS. This fact represents another beneficial outcome of incorporating the cross-derivative con- straints in the process of phase recovery. It should be taken into account that, although the shape of does not change the energy of the PSF i, it plays a crucial role in determining its spatial behaviour In the section that follows, it will be shown that even small inaccuracies in reconstruction of (j) could be translated into dramatic difference in the quality of image deconvolution. B. Image deconvolution As a next step, the phase estimates obtained using the CCS- and DCS-based methods for r — 0.5 were combined with the aperture function A to result in their respective estimates of the PSF i. These estimates were subsequently used to deconvolve (c) (d) Fig. 5. (Subplot (a)) Phase reconstructed obtained by means of CCS for SNR = 40 dB and r = 0.5; (Subplot (b)) Phase reconstructed obtained by means of DCS for the same values of SNR and r.; (Subplot (c) and (d)) Corresponding error maps for CCS and DCS. Fig. 6. MSE of phase reconstruction obtained with different methods as a function of SNR. Here, the dashed and solid lines con'espond to CCS and DCS, respectively, and r = 0.5. a number of test images such as "Satellite", "Saturn", "Moon" and "Galaxy". All the test images were blurred with an original PSF, followed by their contamination with additive Gaussian noise of different levels. As an example. Fig. p\ shows the "Satellite" image (subplot (a)) along with its blurred and noisy version (subplot (b)). Using the PSF estimates, the deconvolution was carried out using the method detailed in 1 1 1 1. For the sake of comparison, the deconvolution was also performed using the PSF recovered from dense sampling (DS) of (p. Note that this reconstruction is expected to have the best accuracy, since it neither involves undersampling nor requires a CS-based phase estimation. All the deconvolved images have been compared with their IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 (a) (b) Fig. 7. Satellite image (subplot (a)) and its blurred and noisy version (subplot (b)). (a) (b) Fig. 8. (Subplot (a)) Image estimate obtained with the CCS-based method for phase recovery (SSIM = 0.917); (Subplot (b)) Image estimate obtained with the DCS-based method for phase recovery (SSIM = 0.781). has been shown to yield phase estimates of substantially better quality as compared to the case of CCS. In this paper, our main focus has been on simplifying the structure of the SHI through reducing the number of its wavefront lenses, while compensating for the effect of undersampling by using the theory of CS augmented by the cross-derivative constraint. The solution was computed using the Bregman algorithm, which provides a computationally ef- ficient framework to carry out the constrained phase recovery. Moreover, the resulting phase estimates were used to recover their associated PSF, which was subsequently used for image deconvolution. It was shown that the DCS-based estimation of (j) with r = 0.3 results in image reconstructions of the quality comparable to that of DS, while substantially outperforming the results obtained with CCS. While the proposed method offers a practical solution to the problem of phase estimation in adaptive optics, some interesting questions about the theoretical aspects of DCS still lay open. In particular, the question of theoretical performance of CS in the presence of side information on the source signal needs to be addressed through future research. Acknowledgment This work was supported in part by the Natural Sciences and Engineering Research Council of Canada and in part by Ontario Early Researcher Award program, which are gratefully acknowledged. The authors would also like to acknowledge Sudipto Dolui for his helpful comments as well as for provid- ing deconvolution codes. original counterparts in terms of PSNR as well as of the structural similarity index (SSIM) of |36|, , which is believed to be a better indicator of perceptual image quality |37|. The resulting values of the comparison metrics are summarized in Table 1, while Fig.lSlshows the deconvolution results produced by the CCS- and DCS-based methods. The above results demonstrate the importance of accurate phase recovery, where even a relatively small phase error can have a dramatic effect on the quality of image deconvolution. Under such conditions, the proposed method produces image reconstructions of a superior quality as compared to the case of CCS. Moreover, comparing the results of Table 1, one can see that DS only slightly outperforms DCS in terms of PSNR and SSIM, while in many practical cases, the difference between the performances of these methods are hard to detect visually. VII. Discussion and conclusions In the present paper, the applicability of DCS to the problem of reconstruction of optical images has been demonstrated. It was shown that, in the presence of atmospheric turbulence, the phase <j) of the GPF P = AeJ'^ is a random function, which needs to be measured using adaptive optics. To simplify the complexity of the latter, a CS-based approach has been proposed. As opposed to CCS, however, the proposed method performs phase reconstruction subject to an additional con- straint, which stems from the property of V0 to be a potential field. The resulting algorithm (referred to as the DCS method) References [1] J. Yang, J. Wright, T. S. Huang, and Y. Ma, "linage super-resolution via sparse representation," IEEE Transactions on Image Processing, vol. 19, pp. 2861-2873, 2010. [2] M. Elad and M. Aharon, "Image denoising via sparse and redundant representations over learned dictionaries," IEEE Transactions on Image Processing, vol. 17. pp. 3736-3745, 2006. [3] IJ. Mairal, G. Sapiro, and M. Elad, "Learning multiscale sparse representations for image and video restoration," Multiscale Modeling and Simulation, vol. 7, pp. 214-241, 2008. [4] R. T. Paul, "Review of robust video watermarking techniques," IJCA Special Issue on Computational Science, , no. 3, pp. 90-95, 2011. [5] G. D. Boreman, Modulation Transfer Function in Optical and Electro- Optical Systems, SPIE Optical Engineering Press, Bellingham, Wash- ington, 2001. [6] R. T. Paul, "Blind deconvolution via cumulant extrema," IEEE Signal Processing Magazine, , no. 3, pp. 2A-A2, 1996. [7] D. Kundur and D. Hatzinakos, "Blind image deconvolution," IEEE Signal Processing Magazine, , no. 3, pp. 43-64, 1996. [8] S. Geman and D. Geman, "Stochastic relaxation, gibbs distribution and the bayesian restoration of images," IEEE Trans. Pattern Analysis and Machine Intelligence, vol. PAMI-6, pp. 721-741, 1984. [9] V. Ton'e T. Poggio and C. Koch, "Computational vision and regulariza- tion theory," Nature, vol. 317, pp. 314-319, 1985. [10] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms," Phys. D, vol. 60, pp. 259-268, November 1992. [11] A. Chambolle, "An algorithm for total variation minimization and applications," Journal of Mathematical Imaging and Vision, vol. 20, pp. 89-97, 2004. [12] O. Michailovich and A. Tannenbaum, "Blind deconvolution of medical ultrasound images: Parametric inverse filtering approach," IEEE Trans- actions on Image Processing, vol. 16, no. 12, pp. 3005-3019, December 2007. [13] W. H. Richardson, "Bayesian-based iterative method of image restora- tion," / Opt. Soc. Am. A, vol. 62, no. 1, pp. 55-59, 1972. IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012 TABLE I SSIM AND PSNR COMPARISONS OF PHASE RECOVERY RESULTS Image Satellite Saturn Moon Galaxy Noise std 10-^ 0.001 0.003 0.005 10--'' 0.001 0.003 0.005 1 10-^ 0.001 0.003 0.005 10-''' 0.001 0.003 0.005 PSNR comparison (in dB) Blurred 14.06 14.06 14.06 14.05 17.78 17.78 17.78 17.78 19.98 19.97 19.97 19.97 18.79 18.79 18.78 18.78 DS 27.97 27.75 25.97 22.43 31.49 31.08 28.50 23.89 25.06 25.04 24.83 21.76 23.58 23.60 23.38 20.93 CS 17.06 16.93 16.54 15.63 23.42 23.38 22.80 20.55 22.36 22.38 22.30 19.73 21.16 21.12 20.64 18.46 DCS 27.42 27.22 25.56 22.22 31.02 30.65 28.30 23.72 25.00 24.99 24.78 21.73 23.52 23.54 23.32 20.86 SSIM comparison Blurred 0.200 0.200 0.199 0.197 0.226 0.226 0.226 0.175 0.512 0.512 0.509 0.504 0.257 0.257 0.257 0.254 DS 0.730 0.720 0.554 0.269 0.688 0.660 0.506 0.228 0.645 0.642 0.607 0.552 0.493 0.495 0.501 0.397 CS 0.349 0.344 0.306 0.206 0.424 0.416 0.348 0.212 0.539 0.538 0.493 0.488 0.348 0.347 0.326 0.224 DCS 0.674 0.667 0.519 0.263 0.656 0.641 0.483 0.223 0.643 0.640 0.604 0.549 0.490 0.491 0.501 0.393 [14] L. B. Lucy, "An iterative technique for the rectification of observed distributions," Astron. /, vol. 79, no. 6, pp. 745-754, 1974. [15] P. A. Jansson, Deconvolution of Images and Spectra, Opt. Eng. 36, 3224, 1997. [16] D. Dayton, B. Pierson, B. Spielbusch, and J. Gonglewski, "Atmospheric structure function measurements with a shack-hartmann wave-front sensor," Journal of Mathematical Imaging and Vision, vol. 20, pp. 89- 97, 2004. [17] R. G. Lane and M. Tallon, "Wave-front reconstruction using a shack- hartmann sensor," Applied Optics, vol. 31, pp. 6902-6908, 1992. [18] M. Hosseini and O. Michailovich, "Derivative compressive sampling with application to phase unwrapping," in Proceedings of EUSIPCO, Glasgow, UK, August 2009. [19] M. C. Roggemann and B. M. Welsh, Imaging Through Turbulence, CRC Press, Boca Raton, FL, 1996. [20] D. L. Fried, "Statistics of a geometric representation of wavefront distortion," J. Opt. Soc. Am., vol. 55, pp. 1427-1431, 1965. [21] O. Michailovich and A. Tannenbaum, "A fast approximation of smooth functions from samples of partial derivatives with application to phase unwrapping," Signal Processing, vol. 88, pp. 358-374, 2008. [22] Y. Tsaig and D. L. Donoho, "Compressed sensing," IEEE Trans. Inform. Theory, vol. 52, pp. 1289-1306, 2006. [23] J. Primot, G. Rousset, and J. C. Fontanella, "Deconvolution from wave- front sensing: a new technique for compensating turbulence-degraded images," / Opt. Soc. Am., vol. 7, pp. 1598-1608, 1990. [24] C. E. Shannon, "Communication in the presence of noise," Proceedings of the IRE, vol. 37, no. 1, pp. 10-21, 1949. [25] E. J. Candes, J. Romberg, and T. Tao, "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency informa- tion," IEEE Trans, hifo. Theory; vol. 52, no. 2, pp. 489-509, 2006. [26] D. L. Donoho and Y. Tsaig, "Fast solution of ii-norm minimization problems when the solution may be sparse," 2006. [27] W. Yin, S. Osher, D. Goldfarb, and J. Darbon, "Bregman iterative al- gorithms for ^ J -minimization with applications to compressed sensing," SIAM J. Imaging Sciences, vol. 1, no. I, pp. 143-168, 2008. [28] J. Cai, S. Osher, and Z. Shen, "Split bregman methods and frame based image restoration," Multiscale Modeling & Simulation, vol. 8, no. 2, pp. 337-369, 2009. [29] S. S. Chen and D. L. Donoho, "Atomic decomposition by basis pursuit," Siam Journal on Scientific Computing, vol. 20, 1998. [30] I. Daubechies, M. Defrise, and C. D. Mol, "An iterative thresholding algorithm for linear inverse problems with a sparsity constraint," Com- munications on Pure and Applied Mathematics, vol. 57, pp. 1413-1457, 2004. [31] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," SIAM Journal on Imaging Sciences, vol. 2, pp. 183-202, 2009. [32] A. Savitzky and M. J. E. Golay, "Smoothing and differentiation of data by simplified least squares procedures," Anal. Chem., vol. 36, pp. 1627-1639, 1964. [33] A. N. Tikhonov and V. Y. Arsenin, "Solutions of ill-posed problems," 1977. [34] J. D. Schmidt, Numerical .simulation of optical wave propagation with examples in MATLAB, SPIE, Washington, 2010. [35] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Reg. Conf. Series in Applied Math. SIAM, 1992. [36] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. R Simoncelli, "Image quality assessment: From error visibility to structural similarity," IEEE Trans. Image Processing, vol. 13, no. 4, pp. 600-612, 2004. [37] Z. Wang and A. C. Bovik, "Mean squared error: love it or leave it? - a new look at signal fidelity measures," IEEE Signal Processing Magazine, vol. 26, no. 1, pp. 98-117, 2009.